program main

    use seakeeping_linalg, only: gemm
    use seakeeping_kinds, only: rk

    interface
        subroutine say_hello() bind(c)
        end subroutine say_hello
    end interface

    integer :: i, j
    real(rk) :: rmat1(2, 2), rmat2(2, 2), rmat3(2, 2)
    complex(rk) :: cmat1(2, 2), cmat2(2, 2), cmat3(2, 2)

    call say_hello()

    rmat1 = 1; rmat2 = 2
    cmat1 = (3, 3); cmat2 = (4, 4); cmat2(2, 2) = (3, 7)

    rmat3 = gemm(rmat1, rmat2); cmat3 = gemm(cmat1, cmat2)

    print '(a)', "Original real matrix:"
    print "(2(4x,es10.3))", ((rmat1(i, j), j=1, 2), i=1, 2)
    print "(2(4x,es10.3))", ((rmat2(i, j), j=1, 2), i=1, 2)
    print '(a)', "Result real matrix:"
    print "(2(4x,es10.3))", ((rmat3(i, j), j=1, 2), i=1, 2)

    print '(a)', "Original complex matrix:"
    print "(4(4x,es10.3))", ((cmat1(i, j), j=1, 2), i=1, 2)
    print "(4(4x,es10.3))", ((cmat2(i, j), j=1, 2), i=1, 2)
    print '(a)', "Result complex matrix:"
    print "(4(4x,es10.3))", ((cmat3(i, j), j=1, 2), i=1, 2)

end program main
